#### FOREIGN FIGHTERS AND CONFLICT-RELATED SEXUAL VIOLENCE. Austin C. Doctor. 2020. International Studies Quarterly. 

library(car)
library(foreign)
library(readstata13)
library(readxl)
library(stargazer)
library(lmtest)
library(clubSandwich)
library(dplyr)
library(plyr) 
library(lubridate)
library(gridExtra)
library(margins)
library(ggplot2)
library(reshape2)
library(MASS)
library(lme4)
library(nnet)

getwd()
setwd("~/Working Papers/Foreign Fighters and CRSV/ISQ Manuscript/Replication Files (Doctor ISQ 2020)")


## UPLOAD DATA ##
data_ff_svac <- read.dta13("replication_doctorISQ_crsvFF.dta")



##################################### MAIN MANUSCRIPT FIGURES #####################################
# Discrete Change in Predicted Probabilities (Figure 1)
plot(0, -0.0893949, ylab= "", xlab="", ylim=c(-0.17, 0.12), xlim=c(-0.3, 1.3), type="p", pch=0, lwd=2, axes=F, main="")
box()
axis(1, at=c(0,0.5,1), labels=c("None", "Some", "High"))
axis(2)
title(ylab="Discrete Change in Predicted Probability", xlab="Prevalence of Sexual Violence")
segments(0, -0.1545031, 0,  -0.0242867, lty=1, col="black")
lines(0.5, 0.0434371, type="p", pch= 0, lwd=2, col="black")
segments(0.5, 0.0050098, 0.5, 0.0818645, lty=1, col="black")
lines(1, 0.0459578, type="p", pch= 0, lwd=2, col="black")
segments(1, 0.0130408, 1, 0.0788748, lty=1, col="black")
abline(h=0, lty=2)



# ISIL Files: Distribution of FF Entrants to ISIL from 2011 to 2014 (Figure 2)
data_isil_files <- read_excel("CGW_FFentries.xlsx")
p <- ggplot(data_isil_files, aes(x=Date, y=NumberFF)) + geom_bar(stat="identity") + coord_cartesian(ylim = c(0,800))
p <- p + labs(title = "", x = "", y="Number of New Foreign Fighter Arrivals") 
p <- p + scale_fill_grey() + theme_classic()
p <- p + scale_x_discrete(breaks=c("2011-01","2012-01","2013-01","2014-01","2014-12"), labels=c("January 2011",  "January 2012", "January 2013",  "January 2014", "December 2014"), expand=c(0,2))
p <- p + theme(axis.title.y = element_text(margin = margin(t = 0, r = 10, b = 0, l = 5)), plot.margin=unit(c(0,1,0,0.5),"cm"))
p



# Distribution of SVAC in Islamic State in Iraq and Levant (Figure 3)
isil_svac <- subset(data_svac_raw, subset=c(actor %in% "IS" & location %in% c("Iraq") & year %in% c(2011:2015)))
isil_svac$svac <- ifelse(isil_svac$state_prev %in% c(2,3), 2, ifelse(isil_svac$ai_prev %in% c(2,3), 2, ifelse(isil_svac$hrw_prev %in% c(2,3), 2, ifelse(isil_svac$state_prev %in% 1, 1, ifelse(isil_svac$ai_prev %in% 1, 1, ifelse(isil_svac$hrw_prev %in% 1, 1, 0))))))

barplot(svac ~ year, data= isil_svac, ylim=c(-0.1,3), ylab= "Prevalence of ISIL Sexual Violence", xlab="Year")

plot(c(2010.5,2015.5), c(-0.3, 2.3), type= "n", xlab = "", ylab = "Prevalence of ISIL Sexual Violence", axes=F)
box()
axis(2,at=seq(0,2,1), labels=c("None", "Some", "High"))
axis(1,at=seq(2011,2015,1))
rect(2010.75, -0.25, 2011.25, 0.25, density = 20, border = "black")
rect(2011.75, -0.25, 2012.25, 0.25, density = 20, border = "black")
rect(2012.75, -0.25, 2013.25, 0.25, density = 20, border = "black")
rect(2013.75, 1.75, 2014.25, 2.25, density = 20, border = "black")
rect(2014.75, 1.75, 2015.25, 2.25, density = 20, border = "black")





##################################### MAIN MANUSCRIPT MODELS #####################################
ologit1 <- polr(as.factor(SVAC_reb_ordinal) ~ foreign_f + time + I(time^2) + I(time^3), data= data_ff_svac, Hess=T); stargazer(ologit1, type="text")
ologit2 <- polr(as.factor(SVAC_reb_ordinal) ~ foreign_f + secessionist + islamist + rel_strength + forced_recruit + loot + time + I(time^2) + I(time^3), data= data_ff_svac, Hess=T); stargazer(ologit2, type="text")
ologit3 <- polr(as.factor(SVAC_reb_ordinal) ~ foreign_f + secessionist + islamist + rel_strength + forced_recruit + loot + duration + govtbestfatal_ln + lnpop + time + I(time^2) + I(time^3), data= data_ff_svac, Hess=T); stargazer(ologit3, type="text")
ologit4 <- polr(as.factor(SVAC_reb_ordinal) ~ foreign_f + SVAC_reb_lag + secessionist + islamist + rel_strength + forced_recruit + loot + time + I(time^2) + I(time^3), data= data_ff_svac, Hess=T); stargazer(ologit4, type="text")
ologit5 <- polr(as.factor(SVAC_reb_ordinal) ~ foreign_f + SVAC_reb_lag + secessionist + islamist + rel_strength + forced_recruit + loot + duration + govtbestfatal_ln + lnpop + time + I(time^2) + I(time^3), data= data_ff_svac, Hess=T); stargazer(ologit5, type="text")

stargazer(ologit1, ologit2, ologit3, ologit4, ologit5, type="text")


## clustered standard errors (conflict dyad)
clval1 <- sandwich::vcovCL(ologit1,factor(data_ff_svac$dyadid)); coeftest(ologit1, vcov= clval1)
clval2 <- sandwich::vcovCL(ologit2,factor(data_ff_svac$dyadid)); coeftest(ologit2, vcov= clval2)
clval3 <- sandwich::vcovCL(ologit3,factor(data_ff_svac$dyadid)); coeftest(ologit3, vcov= clval3)
clval4 <- sandwich::vcovCL(ologit4,factor(data_ff_svac$dyadid)); coeftest(ologit4, vcov= clval4)
clval5 <- sandwich::vcovCL(ologit5,factor(data_ff_svac$dyadid)); coeftest(ologit5, vcov= clval5)







##################################### APPENDIX AND ROBUSTNESS CHECKS #####################################

#### Appendix Table: Descriptive Statistics of Model Variables
stargazer(data_ff_svac, type="text", omit.summary.stat=c("p25", "p75")) # N = 535


#### Appendix Table: Conditional Distribution of Main Variables
table(data_ff_svac$foreign_f, data_ff_svac$SVAC_reb_ordinal)


#### Appendix Figure: Distrubtion of FFs Over Time By Group 
dist.ff <- aggregate(data_ff_svac$dyadid, by=list(data_ff_svac$year, data_ff_svac$foreign_f), FUN = function(x) length(unique(x)))
ff.a <- dist.ff[dist.ff$Group.2 == "0", ]
ff.b <- dist.ff[dist.ff $Group.2 == "1", ]
ff.c <- merge(ff.a, ff.b, by = c("Group.1"), all=T)
ff.c$ratio <- (ff.c$x.y)/(ff.c$x.x + ff.c$x.y); as.data.frame(ff.c); ff.c[is.na(ff.c)] <- 0
ff.c$total <- ff.c$x.x + ff.c$x.y; as.data.frame(ff.c); ff.c[is.na(ff.c)] <- 0
with(ff.c, plot(ratio ~ Group.1, type="l", lwd=1, ylim=c(0,1), main="", ylab="Proportion of Rebel Groups with Foreign Fighters", xlab="Time", lty=1))


#### Appendix Figure: Distrubtion of SVAC Over Time By Group 
dist.svac <- aggregate(data_ff_svac$dyadid, by=list(data_ff_svac$year, data_ff_svac$SVAC_reb), FUN = function(x) length(unique(x)))
svac.a <- dist.svac[dist.svac$Group.2 == "0", ]
svac.b <- dist.svac[dist.svac$Group.2 == "1", ]
svac.c <- merge(svac.a, svac.b, by = c("Group.1"), all=T)
svac.c$ratio <- (svac.c$x.y)/(svac.c$x.x + svac.c$x.y); as.data.frame(svac.c); svac.c[is.na(svac.c)] <- 0
svac.c$total <- svac.c$x.x + svac.c$x.y; as.data.frame(svac.c); svac.c[is.na(svac.c)] <- 0
with(svac.c, plot(ratio ~ Group.1, type="l", lwd=1, ylim=c(0,1), main="", ylab="Proportion of Rebel Groups Perpetrating SVAC", xlab="Time", lty=1))



### Appendix Figure: Distribution of the Dependent Variable
barplot(table(data_ff_svac$SVAC_reb_ordinal), main="Distribution of Conflict-Related Sexual Violence", ylim=c(0,500)); table(data_ff_svac$SVAC_reb_ordinal)
table(data_ff_svac$SVAC_reb_ordinal)


#### Appendix Table: Original 4-Point Measure of SVAC ####
data_ff_svac$SVAC_reb_OG <- ifelse(data_ff_svac$state_prev %in% c(3), 3, ifelse(data_ff_svac$ai_prev %in% c(3), 3, ifelse(data_ff_svac$hrw_prev %in% c(3), 3, ifelse(data_ff_svac$state_prev %in% c(2), 2, ifelse(data_ff_svac$ai_prev %in% c(2), 2, ifelse(data_ff_svac$hrw_prev %in% c(2), 2, ifelse(data_ff_svac$state_prev %in% 1, 1, ifelse(data_ff_svac$ai_prev %in% 1, 1, ifelse(data_ff_svac$hrw_prev %in% 1, 1, 0)))))))))

table(data_ff_svac$SVAC_reb_OG)

ologit1_OG <- polr(as.factor(SVAC_reb_OG) ~ foreign_f + time + I(time^2) + I(time^3), data= data_ff_svac, Hess=T)
ologit2_OG <- polr(as.factor(SVAC_reb_OG) ~ foreign_f + secessionist + islamist + rel_strength + forced_recruit + loot + SVAC_reb_lag + time + I(time^2) + I(time^3), data= data_ff_svac, Hess=T)
ologit3_OG <- polr(as.factor(SVAC_reb_OG) ~ foreign_f + secessionist + islamist + rel_strength + forced_recruit + loot + duration + govtbestfatal_ln + lnpop + SVAC_reb_lag + time + I(time^2) + I(time^3), data= data_ff_svac, Hess=T)

stargazer(ologit1_OG, ologit2_OG, ologit3_OG, type="text")



#### Appendix Table: Additional Control Variables
data_reb_gov <- read.dta13("replication_doctorISQ_data_rebgov.dta")
stargazer(data_reb_gov, type="text")

ologit1_add <- polr(as.factor(SVAC_reb_ordinal) ~ foreign_f + SVAC_gov_binary + SVAC_reb_lag + time + I(time^2) + I(time^3), data= data_reb_gov, Hess=T)
ologit2_add <- polr(as.factor(SVAC_reb_ordinal) ~ foreign_f + secessionist + islamist + rel_strength + forced_recruit + + loot + rebelsupport + terrcont + centcontrol + SVAC_gov_binary + SVAC_reb_lag + time + I(time^2) + I(time^3), data= data_reb_gov, Hess=T)
ologit3_add <- polr(as.factor(SVAC_reb_ordinal) ~ foreign_f + secessionist + islamist + rel_strength + forced_recruit + loot + rebelsupport + terrcont + centcontrol + SVAC_gov_binary + intensity + duration + govtbestfatal_ln + lnpop + SVAC_reb_lag + time + I(time^2) + I(time^3), data= data_reb_gov, Hess=T)

stargazer(ologit1_add, ologit2_add, ologit3_add, type="text")



#### Appendix Table: Country-Level Fixed Effects ####
fe_ologit1 <- polr(as.factor(SVAC_reb_ordinal) ~ foreign_f + time + I(time^2) + I(time^3) + as.factor(location), data= data_ff_svac, Hess=T)
fe_ologit2 <- polr(as.factor(SVAC_reb_ordinal) ~ foreign_f + secessionist + islamist + rel_strength + forced_recruit + loot + SVAC_reb_lag + time + I(time^2) + I(time^3) + as.factor(location), data= data_ff_svac, Hess=T)
fe_ologit3 <- polr(as.factor(SVAC_reb_ordinal) ~ foreign_f + secessionist + islamist + rel_strength + forced_recruit + loot + duration + govtbestfatal_ln + lnpop + SVAC_reb_lag + time + I(time^2) + I(time^3) + as.factor(location), data= data_ff_svac, Hess=T)

stargazer(fe_ologit1, fe_ologit2, fe_ologit3, type="text")


#### Appendix Table: Estimate All Models on Sample Aggregated to the Conflict-Dyad Level
#See STATA do file

